clear all
set maxvar 10000
set matsize 11000
set more off


use Setting_Sun_dataset.dta, clear

sort hour
gen id = _n
tsset id
tab month, gen(m)

*******************************************************************************************************
** Average change in hourly RTM price (January 1, 2013 through May 31, 2017) -- Equation (1).
** Note: the coefficients on daily_solar, daily_wind, and spot are displayed in Figure 2. The
** coefficients on daily_solar and daily_wind also appear in Figures A4, A9, and A14.


forvalues h = 1/24 {
	quietly: newey rtm_total_lmp daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..3]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..3, 1..3]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename beta_lmp3 beta_spot
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind
rename stddev_lmp3 stddev_spot

export excel hour-stddev_spot using "Regression_Output.xlsx", sheet("RTM_Price_Change") firstrow(varlabels) cell(A1) sheetmodify

*******************************************************************************************************	
** Average change in hourly quantity supplied to CAISO by source.
** Note: the coefficents on daily_solar are displayed in the top panel of Figure 3. The coefficients on
** daily_wind are displayed in Figure A6.


foreach i in thermal large_hydro nuclear imports other {
use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)
gen other = renewables - solar - wind
forvalues h = 1/24 {
	quietly: newey `i' daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_`h' = beta_raw[1,1..2]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..2, 1..2]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_`h' = vecdiag(std_raw)
}

matrix define beta = (beta_1)
matrix define stddev = (stddev_1)

forvalues h = 2/24 {
	matrix define beta = (beta \ beta_`h')
	matrix define stddev = (stddev \ stddev_`h')
}

clear
set obs 24
gen hour = _n
svmat beta
svmat stddev
rename beta1 beta_solar
rename beta2 beta_wind
rename stddev1 stddev_solar
rename stddev2 stddev_wind

export excel hour-stddev_wind using "Regression_Output.xlsx", sheet("`i'_Change") firstrow(varlabels) cell(A1) sheetmodify
}


*******************************************************************************************************	
** Average change in net generation from CAISO fossil-fuel units in the EPA's CEMS data.
** Note: the coefficents on daily_solar are displayed in the bottom panel of Figure 3 and Figure A17.


foreach i in ccgt gt st {
use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)
forvalues h = 1/24 {
	quietly: newey net_load_`i'_caiso daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_`h' = beta_raw[1,1..2]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..2, 1..2]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_`h' = vecdiag(std_raw)
}

matrix define beta = (beta_1)
matrix define stddev = (stddev_1)

forvalues h = 2/24 {
	matrix define beta = (beta \ beta_`h')
	matrix define stddev = (stddev \ stddev_`h')
}

clear
set obs 24
gen hour = _n
svmat beta
svmat stddev
rename beta1 beta_solar
rename beta2 beta_wind
rename stddev1 stddev_solar
rename stddev2 stddev_wind

export excel hour beta_solar stddev_solar using "Regression_Output.xlsx", sheet("`i'_Change") firstrow(varlabels) cell(A1) sheetmodify
}


*******************************************************************************************************	
** Average change in daily CO2 and NOx emisions by season.
** Note: the estimates of the average daily change in CAISO emissions, large hydro output, and imports
** are displayed in Table 3. The coefficients of the average change in hourly CAISO CO2 emissions from the model
** below are displayed in Figure 4. The coefficients on the average change in hourly large-hydro output and 
** net-imports into CAISO by season are displayed in Figure 5.

use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)

gen season = 1
replace season = 2 if month >=3 & month <=5
replace season = 3 if month >=6 & month <=8
replace season = 4 if month >=9 & month <=11
tab hour, gen(h)
save "seasonal_dataset_temp.dta", replace

forvalues s = 1/4 {

use "seasonal_dataset_temp.dta", clear
keep if season == `s'

forvalues h = 1/24 {
	gen daily_solar_h`h' = daily_solar * h`h'
}
forvalues h = 1/24 {
	gen daily_wind_h`h' = daily_wind * h`h'
}
forvalues h = 1/24 {
	gen spot_h`h' = spot * h`h'
}
forvalues h = 1/24 {
	gen inches_lag_h`h' = inches_lag * h`h'
}
forvalues h = 1/24 {
	gen load_h`h' = load * h`h'
}
forvalues h = 1/24 {
	forvalues m = 1/12 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

gen hydro_imports = large_hydro + imports

sort year month day hour
gen id2 = _n
tsset id2

quietly: newey co2_caiso daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24, lag(168) force

lincom daily_solar_h1 + daily_solar_h2 + daily_solar_h3 + daily_solar_h4 + daily_solar_h5 + daily_solar_h6 + daily_solar_h7 + daily_solar_h8 + daily_solar_h9 + daily_solar_h10 + daily_solar_h11 + daily_solar_h12 + daily_solar_h13 + daily_solar_h14 + daily_solar_h15 + daily_solar_h16 + daily_solar_h17 + daily_solar_h18 + daily_solar_h19 + daily_solar_h20 + daily_solar_h21 + daily_solar_h22 + daily_solar_h23 + daily_solar_h24
scalar define co2_s`s' = r(estimate)
scalar define co2_s`s'_std = r(se)


quietly: newey large_hydro daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24, lag(168) force

lincom daily_solar_h1 + daily_solar_h2 + daily_solar_h3 + daily_solar_h4 + daily_solar_h5 + daily_solar_h6 + daily_solar_h7 + daily_solar_h8 + daily_solar_h9 + daily_solar_h10 + daily_solar_h11 + daily_solar_h12 + daily_solar_h13 + daily_solar_h14 + daily_solar_h15 + daily_solar_h16 + daily_solar_h17 + daily_solar_h18 + daily_solar_h19 + daily_solar_h20 + daily_solar_h21 + daily_solar_h22 + daily_solar_h23 + daily_solar_h24
scalar define hydro_s`s' = r(estimate)
scalar define hydro_s`s'_std = r(se)


quietly: newey imports daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24, lag(168) force

lincom daily_solar_h1 + daily_solar_h2 + daily_solar_h3 + daily_solar_h4 + daily_solar_h5 + daily_solar_h6 + daily_solar_h7 + daily_solar_h8 + daily_solar_h9 + daily_solar_h10 + daily_solar_h11 + daily_solar_h12 + daily_solar_h13 + daily_solar_h14 + daily_solar_h15 + daily_solar_h16 + daily_solar_h17 + daily_solar_h18 + daily_solar_h19 + daily_solar_h20 + daily_solar_h21 + daily_solar_h22 + daily_solar_h23 + daily_solar_h24
scalar define imports_s`s' = r(estimate)
scalar define imports_s`s'_std = r(se)


quietly: newey hydro_imports daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24, lag(168) force

lincom daily_solar_h1 + daily_solar_h2 + daily_solar_h3 + daily_solar_h4 + daily_solar_h5 + daily_solar_h6 + daily_solar_h7 + daily_solar_h8 + daily_solar_h9 + daily_solar_h10 + daily_solar_h11 + daily_solar_h12 + daily_solar_h13 + daily_solar_h14 + daily_solar_h15 + daily_solar_h16 + daily_solar_h17 + daily_solar_h18 + daily_solar_h19 + daily_solar_h20 + daily_solar_h21 + daily_solar_h22 + daily_solar_h23 + daily_solar_h24
scalar define hydro_imports_s`s' = r(estimate)
scalar define hydro_imports_s`s'_std = r(se)


quietly: newey nox_caiso daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24, lag(168) force

lincom daily_solar_h1 + daily_solar_h2 + daily_solar_h3 + daily_solar_h4 + daily_solar_h5 + daily_solar_h6 + daily_solar_h7 + daily_solar_h8 + daily_solar_h9 + daily_solar_h10 + daily_solar_h11 + daily_solar_h12 + daily_solar_h13 + daily_solar_h14 + daily_solar_h15 + daily_solar_h16 + daily_solar_h17 + daily_solar_h18 + daily_solar_h19 + daily_solar_h20 + daily_solar_h21 + daily_solar_h22 + daily_solar_h23 + daily_solar_h24
scalar define nox_s`s' = r(estimate)
scalar define nox_s`s'_std = r(se)


matrix define daily_impact_s`s' = [co2_s`s', hydro_s`s', imports_s`s', hydro_imports_s`s', nox_s`s']
matrix define daily_impact_std_s`s' = [co2_s`s'_std, hydro_s`s'_std, imports_s`s'_std, hydro_imports_s`s'_std, nox_s`s'_std]

}

** Note: the daily impact values represent the average change in CAISO CO2, NOx, large-hydro, and net imports per GWh of daily solar.
** These are values are displayed in Table 3.

matrix define daily_impact = [daily_impact_s1 \ daily_impact_s2 \ daily_impact_s3 \ daily_impact_s4]
matrix define daily_impact_std = [daily_impact_std_s1 \ daily_impact_std_s2 \ daily_impact_std_s3 \ daily_impact_std_s4]


clear
svmat daily_impact
svmat daily_impact_std

** Rows 1 through 4 correspond to the estimates from seasons 1 through 4.

rename daily_impact1 co2
rename daily_impact2 hydro
rename daily_impact3 imports
rename daily_impact4 hydro_imports
rename daily_impact5 nox
rename daily_impact_std1 co2_std
rename daily_impact_std2 hydro_std
rename daily_impact_std3 imports_std
rename daily_impact_std4 hydro_imports_std
rename daily_impact_std5 nox_std

export excel co2-nox_std using "Regression_Output.xlsx", sheet("Emission_Impacts") firstrow(varlabels) cell(A1) sheetmodify



*******************************************************************************************************	
** The following section predicts counterfactual 2016 RTM prices for different levels of installed solar
** capacity. Using the predicted counterfactual prices, we then predict the operating profits for
** hypothetical conventional generators with a range of different marginal generation costs.

use "seasonal_dataset_temp.dta", clear

** Create daily wind and solar potential capacity factor measures
gen solar_cap_factor = solar_potential / (24 * solar_cap/1000)
gen wind_cap_factor = wind_potential / (24 * wind_cap/1000)

** Create hourly interactions for subsequent DAM price model
forvalues h = 1/24 {
	gen daily_solar_h`h' = 0
	replace daily_solar_h`h' = daily_solar if hour == `h'
}

forvalues h = 1/24 {
	gen daily_wind_h`h' = 0
	replace daily_wind_h`h' = daily_wind if hour == `h'
}

forvalues h = 1/24 {
	gen load_h`h' = 0
	replace load_h`h' = load if hour == `h'
}

forvalues h = 1/24 {
	gen spot_h`h' = 0
	replace spot_h`h' = spot if hour == `h'
}

forvalues h = 1/24 {
	gen inches_lag_h`h' = 0
	replace inches_lag_h`h' = inches_lag if hour == `h'
}

forvalues m = 1/12 {
	forvalues h = 1/24 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

** Create countefactual 2016 datasets with 2, 6, & 10 GWs of installed solar capacity

save "full_temp.dta", replace
gen solar_sample = "actual"

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 2000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_2k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 6000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_6k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 10000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_10k" if solar_sample == ""

** Estimate Curtailed Solar as a function of the daily load, solar potential, wind, etc.
** Note: Results from tobit model are presented in Table A3

gen solar_load = solar_potential * daily_load
tobit solar_curtailed solar_potential solar_load daily_wind spot daily_load inches_lag m1-m11 if hour == 12 & solar_sample == "actual", ll(0)
predict solar_curtailed_fitted, ystar(0,.)

** Counterfactual daily solar is equal to the potential output minus predicted curtailment
replace daily_solar = solar_potential - solar_curtailed_fitted if solar_sample != "actual"

forvalues h = 1/24 {
	replace daily_solar_h`h' = daily_solar if hour == `h' & solar_sample != "actual"
}

** Predict counterfactual RTM prices
** Note: coefficients of average hourly impact of daily solar & wind by season are displayed in Figures A15 & A16.

gen rtm_fitted = .
gen resid = .

forvalues s = 1/4 {
	quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 h1-h24 if solar_sample == "actual" & season == `s'
	predict rtm_fitted_`s'
	predict resid_`s', resid
	replace rtm_fitted = rtm_fitted_`s' if season == `s'
	replace resid = resid_`s' if solar_sample == "actual" & season == `s'
}

replace resid = resid[_n-38685] if solar_sample == "solar_2k"
replace resid = resid[_n-38685] if solar_sample == "solar_6k"
replace resid = resid[_n-38685] if solar_sample == "solar_10k"

gen rtm_hat = rtm_fitted + resid
drop if year != 2016
save "counterfactual_price_by_solar_temp.dta", replace

** Generate 2016 average hourly price profile by solar capacity
** The avg. counterfactual prices are displayed in the top panel of Figures 6 and A7.

reg rtm_hat h1-h24 if solar_sample == "actual", noc
matrix define beta_actual = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_2k", noc
matrix define beta_2k = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_6k", noc
matrix define beta_6k = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_10k", noc
matrix define beta_10k = e(b)'

matrix define avg_price = (beta_actual, beta_2k, beta_6k, beta_10k)
clear
set obs 24
gen hour = _n
svmat avg_price
rename avg_price1 actual
rename avg_price2 solar_2k
rename avg_price3 solar_6k
rename avg_price4 solar_10k

export excel hour-solar_10k using "Counterfactual.xlsx", sheet("Price_Solar") firstrow(varlabels) sheetmodify cell(A1)

** Generate kernel-smoothed distribution of predicted hourly prices
** Note: distributions are displayed in Figure A18

use "counterfactual_price_by_solar_temp.dta", clear
gen p = _n - 6
replace p = . if p > 110

quietly: kdensity rtm_hat if solar_sample == "actual", at(p) gen(dist_actual)
quietly: kdensity rtm_hat if solar_sample == "solar_2k", at(p) gen(dist_2k)
quietly: kdensity rtm_hat if solar_sample == "solar_6k", at(p) gen(dist_6k)
quietly: kdensity rtm_hat if solar_sample == "solar_10k", at(p) gen(dist_10k)

keep p dist_actual dist_2k dist_6k dist_10k
drop if p == .

export excel p-dist_10k using "Counterfactual.xlsx", sheet("Distribution_Solar") firstrow(varlabels) sheetmodify cell(A1)

** Predict share of hours operating and operating profits for hypothetical technologies with different marginal costs.
** Note: the predicted changes in operating profits are displayed in the top panel of Figure 7. The changes in the
** predicted operation time are displayed in the top panel of Figure A19.

use "counterfactual_price_by_solar_temp.dta", clear

forvalues i = 0/100 {

	gen mc_`i' = 0
	replace mc_`i' = 1 if rtm_hat >= `i'
	gen pi_`i' = 0
	replace pi_`i' = rtm_hat - `i' if mc_`i' == 1
	
}

order mc_0-pi_100, sequential
collapse (sum) mc_0-pi_100, by(solar_sample)

mkmat mc_0-mc_100, matrix(operate)
mkmat pi_0-pi_100, matrix(profit)
matrix define profit = profit'
matrix define operate = operate'

clear
set obs 101
gen mc = _n - 1
svmat operate
svmat profit

foreach i in operate profit {
	rename `i'1 `i'_actual
	rename `i'2 `i'_10k
	rename `i'3 `i'_2k
	rename `i'4 `i'_6k
}

export excel mc-profit_6k using "Counterfactual.xlsx", sheet("Profit_Solar") firstrow(varlabels) sheetmodify cell(A1)


*******************************************************************************************************	
** The following section predicts counterfactual 2016 RTM prices for different levels of installed wind
** capacity. Using the predicted counterfactual prices, we then predict the operating profits for
** hypothetical conventional generators with a range of different marginal generation costs.


use "seasonal_dataset_temp.dta", clear

** Create daily wind and solar potential capacity factor measures
gen solar_cap_factor = solar_potential / (24 * solar_cap/1000)
gen wind_cap_factor = wind_potential / (24 * wind_cap/1000)

** Create hourly interactions for subsequent DAM price model
forvalues h = 1/24 {
	gen daily_solar_h`h' = 0
	replace daily_solar_h`h' = daily_solar if hour == `h'
}

forvalues h = 1/24 {
	gen daily_wind_h`h' = 0
	replace daily_wind_h`h' = daily_wind if hour == `h'
}

forvalues h = 1/24 {
	gen load_h`h' = 0
	replace load_h`h' = load if hour == `h'
}

forvalues h = 1/24 {
	gen spot_h`h' = 0
	replace spot_h`h' = spot if hour == `h'
}

forvalues h = 1/24 {
	gen inches_lag_h`h' = 0
	replace inches_lag_h`h' = inches_lag if hour == `h'
}

forvalues m = 1/12 {
	forvalues h = 1/24 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

** Create countefactual 2016 datasets with 1k, 3k, and 6k MW of installed wind capacity

save "full_temp.dta", replace
gen wind_sample = "actual"

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 1000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_1k" if wind_sample == ""

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 3000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_3k" if wind_sample == ""

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 6000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_6k" if wind_sample == ""

** Estimate Curtailed wind as a function of the daily load, wind potential, solar etc.

gen wind_load = wind_potential * daily_load
tobit wind_curtailed daily_solar wind_load wind_potential spot daily_load inches_lag m1-m11 if hour == 12 & wind_sample == "actual", ll(0)
predict wind_curtailed_fitted, ystar(0,.)

** Counterfactual daily wind is equal to the potential output minus predicted curtailment
replace daily_wind = wind_potential - wind_curtailed_fitted if wind_sample != "actual"

forvalues h = 1/24 {
	replace daily_wind_h`h' = daily_wind if hour == `h' & wind_sample != "actual"
}

** Predict counterfactual RTM prices

gen rtm_fitted = .
gen resid = .

forvalues s = 1/4 {
	quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 h1-h24 if wind_sample == "actual" & season == `s'
	predict rtm_fitted_`s'
	predict resid_`s', resid
	replace rtm_fitted = rtm_fitted_`s' if season == `s'
	replace resid = resid_`s' if wind_sample == "actual" & season == `s'
}

replace resid = resid[_n-38685] if wind_sample == "wind_1k"
replace resid = resid[_n-38685] if wind_sample == "wind_3k"
replace resid = resid[_n-38685] if wind_sample == "wind_6k"

gen rtm_hat = rtm_fitted + resid
drop if year != 2016
save "counterfactual_price_by_wind_temp.dta", replace

** Generate 2016 average hourly price profile by wind capacity
** Note: the average hourly prices are displayed in the bottom panels of Figures 6 and A7.

reg rtm_hat h1-h24 if wind_sample == "actual", noc
matrix define beta_actual = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_1k", noc
matrix define beta_1k = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_3k", noc
matrix define beta_3k = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_6k", noc
matrix define beta_6k = e(b)'

matrix define avg_price = (beta_actual, beta_1k, beta_3k, beta_6k)
clear
set obs 24
gen hour = _n
svmat avg_price
rename avg_price1 actual
rename avg_price2 wind_1k
rename avg_price3 wind_3k
rename avg_price4 wind_6k

export excel hour-wind_6k using "Counterfactual.xlsx", sheet("Price_Wind") firstrow(varlabels)  sheetmodify cell (A1)

** Predict share of hours operating and operating profits for hypothetical technologies with different marginal costs.
** Note: the predicted changes in operating profits are displayed in the bottom panel of Figure 7. The changes in the
** predicted operation time are displayed in the bottom panel of Figure A19.

use "counterfactual_price_by_wind_temp.dta", clear

forvalues i = 0/100 {

	gen mc_`i' = 0
	replace mc_`i' = 1 if rtm_hat >= `i'
	gen pi_`i' = 0
	replace pi_`i' = rtm_hat - `i' if mc_`i' == 1
	
}

order mc_0-pi_100, sequential
collapse (sum) mc_0-pi_100, by(wind_sample)

mkmat mc_0-mc_100, matrix(operate)
mkmat pi_0-pi_100, matrix(profit)
matrix define profit = profit'
matrix define operate = operate'

clear
set obs 101
gen mc = _n - 1
svmat operate
svmat profit

foreach i in operate profit {
	rename `i'1 `i'_actual
	rename `i'2 `i'_1k
	rename `i'3 `i'_3k
	rename `i'4 `i'_6k
}

export excel mc-profit_6k using "Counterfactual.xlsx", sheet("Profit_Wind") firstrow(varlabels) sheetmodify cell (A1)



*******************************************************************************************************	
** The following section predicts counterfactual 2016 RTM prices for different levels of installed solar
** capacity imposing the assumption that curtailment levels do not vary with renewable capacity.
** The comparison of the counterfactual prices with and without endogenous curtailment is displayed in
** in the top panel of Figure A8.

use "seasonal_dataset_temp.dta", clear

** Create daily wind and solar potential capacity factor measures
gen solar_cap_factor = solar_potential / (24 * solar_cap/1000)
gen wind_cap_factor = wind_potential / (24 * wind_cap/1000)

** Create hourly interactions for subsequent DAM price model
forvalues h = 1/24 {
	gen daily_solar_h`h' = 0
	replace daily_solar_h`h' = daily_solar if hour == `h'
}

forvalues h = 1/24 {
	gen daily_wind_h`h' = 0
	replace daily_wind_h`h' = daily_wind if hour == `h'
}

forvalues h = 1/24 {
	gen load_h`h' = 0
	replace load_h`h' = load if hour == `h'
}

forvalues h = 1/24 {
	gen spot_h`h' = 0
	replace spot_h`h' = spot if hour == `h'
}

forvalues h = 1/24 {
	gen inches_lag_h`h' = 0
	replace inches_lag_h`h' = inches_lag if hour == `h'
}

forvalues m = 1/12 {
	forvalues h = 1/24 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

** Create countefactual 2016 datasets with 2, 6, & 10 GWs of installed solar capacity

save "full_temp.dta", replace
gen solar_sample = "actual"

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 2000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_2k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 6000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_6k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 10000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_10k" if solar_sample == ""

** Counterfactual daily solar is equal to the potential output minus ACTUAL curtailment
replace daily_solar = solar_potential - solar_curtailed if solar_sample != "actual"

forvalues h = 1/24 {
	replace daily_solar_h`h' = daily_solar if hour == `h' & solar_sample != "actual"
}

** Predict counterfactual RTM prices

gen rtm_fitted = .
gen resid = .

forvalues s = 1/4 {
	quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 h1-h24 if solar_sample == "actual" & season == `s'
	predict rtm_fitted_`s'
	predict resid_`s', resid
	replace rtm_fitted = rtm_fitted_`s' if season == `s'
	replace resid = resid_`s' if solar_sample == "actual" & season == `s'
}

replace resid = resid[_n-38685] if solar_sample == "solar_2k"
replace resid = resid[_n-38685] if solar_sample == "solar_6k"
replace resid = resid[_n-38685] if solar_sample == "solar_10k"

gen rtm_hat = rtm_fitted + resid
drop if year != 2016

** Generate 2016 average hourly price profile by solar capacity

reg rtm_hat h1-h24 if solar_sample == "actual", noc
matrix define beta_actual = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_2k", noc
matrix define beta_2k = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_6k", noc
matrix define beta_6k = e(b)'
reg rtm_hat h1-h24 if solar_sample == "solar_10k", noc
matrix define beta_10k = e(b)'

matrix define avg_price = (beta_actual, beta_2k, beta_6k, beta_10k)
clear
set obs 24
gen hour = _n
svmat avg_price
rename avg_price1 actual
rename avg_price2 solar_2k
rename avg_price3 solar_6k
rename avg_price4 solar_10k

export excel hour-solar_10k using "Counterfactual.xlsx", sheet("Price_Solar_Exogenous") firstrow(varlabels) sheetmodify cell(A1)




*******************************************************************************************************	
** The following section predicts counterfactual 2016 RTM prices for different levels of installed wind
** capacity imposing the assumption that curtailment levels do not vary with renewable capacity.
** The comparison of the counterfactual prices with and without endogenous curtailment is displayed in
** in the bottom panel of Figure A8.

use "seasonal_dataset_temp.dta", clear

** Create daily wind and solar potential capacity factor measures
gen solar_cap_factor = solar_potential / (24 * solar_cap/1000)
gen wind_cap_factor = wind_potential / (24 * wind_cap/1000)

** Create hourly interactions for subsequent DAM price model
forvalues h = 1/24 {
	gen daily_solar_h`h' = 0
	replace daily_solar_h`h' = daily_solar if hour == `h'
}

forvalues h = 1/24 {
	gen daily_wind_h`h' = 0
	replace daily_wind_h`h' = daily_wind if hour == `h'
}

forvalues h = 1/24 {
	gen load_h`h' = 0
	replace load_h`h' = load if hour == `h'
}

forvalues h = 1/24 {
	gen spot_h`h' = 0
	replace spot_h`h' = spot if hour == `h'
}

forvalues h = 1/24 {
	gen inches_lag_h`h' = 0
	replace inches_lag_h`h' = inches_lag if hour == `h'
}

forvalues m = 1/12 {
	forvalues h = 1/24 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

** Create countefactual 2016 datasets with 1k, 3k, and 6k MW of installed wind capacity

save "full_temp.dta", replace
gen wind_sample = "actual"

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 1000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_1k" if wind_sample == ""

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 3000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_3k" if wind_sample == ""

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * 6000 * (24/1000) if wind_sample == ""
replace wind_sample = "wind_6k" if wind_sample == ""

** Counterfactual daily wind is equal to the potential output minus ACTUAL curtailment
replace daily_wind = wind_potential - wind_curtailed if wind_sample != "actual"

forvalues h = 1/24 {
	replace daily_wind_h`h' = daily_wind if hour == `h' & wind_sample != "actual"
}

** Predict counterfactual RTM prices

gen rtm_fitted = .
gen resid = .

forvalues s = 1/4 {
	quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 h1-h24 if wind_sample == "actual" & season == `s'
	predict rtm_fitted_`s'
	predict resid_`s', resid
	replace rtm_fitted = rtm_fitted_`s' if season == `s'
	replace resid = resid_`s' if wind_sample == "actual" & season == `s'
}

replace resid = resid[_n-38685] if wind_sample == "wind_1k"
replace resid = resid[_n-38685] if wind_sample == "wind_3k"
replace resid = resid[_n-38685] if wind_sample == "wind_6k"

gen rtm_hat = rtm_fitted + resid
drop if year != 2016

** Generate 2016 average hourly price profile by wind capacity

reg rtm_hat h1-h24 if wind_sample == "actual", noc
matrix define beta_actual = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_1k", noc
matrix define beta_1k = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_3k", noc
matrix define beta_3k = e(b)'
reg rtm_hat h1-h24 if wind_sample == "wind_6k", noc
matrix define beta_6k = e(b)'

matrix define avg_price = (beta_actual, beta_1k, beta_3k, beta_6k)
clear
set obs 24
gen hour = _n
svmat avg_price
rename avg_price1 actual
rename avg_price2 wind_1k
rename avg_price3 wind_3k
rename avg_price4 wind_6k

export excel hour-wind_6k using "Counterfactual.xlsx", sheet("Price_Wind_Exogenous") firstrow(varlabels)  sheetmodify cell (A1)


*******************************************************************************************************	
** The following section predicts the marginal energy value (during 2016) from an additional MW of
** solar capacity. The predicted changes in the market value of the additional solar output are
** displayed in the top panel of Figure 8.

use "full_temp.dta", clear

** Create hourly capacity factors (net of curtailments)
gen hourly_solar_cap_factor = solar/solar_cap
gen hourly_wind_cap_factor = wind/wind_cap

** Create countefactual 2016 datasets with 2k up to 10k MW of installed solar capacity (500 MW steps)
save "full_temp.dta", replace
gen solar_sample = "actual"
gen sample_id = 0

forvalues i = 2000(500)10000 {
	append using "full_temp.dta"
	replace solar_potential = solar_cap_factor * `i' * (24/1000) if solar_sample == ""
	replace sample_id = `i' if solar_sample == ""
	replace solar_sample = "solar_`i'" if solar_sample == ""
}

** Estimate Curtailed Solar as a function of the daily load, solar potential, wind, etc.

gen solar_load = solar_potential * daily_load
tobit solar_curtailed solar_potential solar_load daily_wind spot daily_load inches_lag m1-m11 if hour == 12 & solar_sample == "actual", ll(0)
predict solar_curtailed_fitted, ystar(0,.)

** Counterfactual daily solar is equal to the potential output minus predicted curtailment
replace daily_solar = solar_potential - solar_curtailed_fitted if solar_sample != "actual"
forvalues h = 1/24 {
	replace daily_solar_h`h' = daily_solar if hour == `h' & solar_sample != "actual"
}

** Predict counterfactual RTM prices

quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24 if solar_sample == "actual"
predict rtm_fitted
predict resid, resid

forvalues i = 2000(500)10000 {
	replace resid = resid[_n-38685] if solar_sample == "solar_`i'"
}

gen rtm_hat = rtm_fitted + resid
drop if year != 2016

** Calculate 2016 revenue from marginal solar capacity

gen marginal_revenue = rtm_hat * hourly_solar_cap_factor
collapse (sum) marginal_revenue, by(sample_id)
drop if sample_id == 0
rename sample_id solar_capacity

export excel solar_capacity marginal_revenue using "Counterfactual.xlsx", sheet("Solar_MR") firstrow(varlabels) sheetmodify cell(A1)


*******************************************************************************************************	
** The following section predicts the marginal energy value (during 2016) from an additional MW of
** wind capacity. The predicted changes in the market value of the additional solar output are
** displayed in the bottom panel of Figure 8.

** Create countefactual 2016 datasets with 1k up to 6k MW of installed wind capacity (500 MW steps)
use "full_temp.dta", clear
gen wind_sample = "actual"
gen sample_id = 0

forvalues i = 1000(500)6000 {

append using "full_temp.dta"
replace wind_potential = wind_cap_factor * `i' * (24/1000) if wind_sample == ""
replace sample_id = `i' if wind_sample == ""
replace wind_sample = "wind_`i'" if wind_sample == ""

}

** Estimate Curtailed Wind as a function of the daily load, wind potential, solar, etc.

gen wind_load = wind_potential * daily_load
tobit wind_curtailed daily_solar wind_load wind_potential spot daily_load inches_lag m1-m11 if hour == 12 & wind_sample == "actual", ll(0)
predict wind_curtailed_fitted, ystar(0,.)

** Counterfactual daily wind is equal to the potential output minus predicted curtailment
replace daily_wind = wind_potential - wind_curtailed_fitted if wind_sample != "actual"
forvalues h = 1/24 {
	replace daily_wind_h`h' = daily_wind if hour == `h' & wind_sample != "actual"
}

** Predict counterfactual RTM prices

quietly: reg rtm_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 m1_h1-m12_h24 if wind_sample == "actual"
predict rtm_fitted
predict resid, resid

forvalues i = 1000(500)6000 {
	replace resid = resid[_n-38685] if wind_sample == "wind_`i'"
}

gen rtm_hat = rtm_fitted + resid
drop if year != 2016

** Calculate 2016 revenue from marginal wind capacity

gen marginal_revenue = rtm_hat * hourly_wind_cap_factor
collapse (sum) marginal_revenue, by(sample_id)
drop if sample_id == 0
rename sample_id wind_capacity

export excel wind_capacity marginal_revenue using "Counterfactual.xlsx", sheet("Wind_MR") firstrow(varlabels) sheetmodify cell(A1)


*****************************************************************************************
*****************************************************************************************
** ADDITIONAL ROBUSTNESS CHECKS

*****************************************************************************************
** The section below estimates the average change in hourly RTM price by using the potential daily renewable output
** as instruments for daily solar and wind generation. For comparison, the standard OLS estimates of Eq. 1 are also
** repeated below over the same time period for which the potential wind and solar generaiton is observed. The
** resulting IV and OLS point estimates are displayed in Figure A12.

use Setting_Sun_dataset.dta, clear
tab month, gen(m)
gen week = week(julian)
gen year_week = year*100 + week
forvalues h = 1/24 {
	quietly: ivreg2 rtm_total_lmp spot load inches_lag m1-m11 (daily_solar daily_wind = solar_potential wind_potential) if hour == `h', robust cluster(year_week)
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..2]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..2, 1..2]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind

export excel hour-stddev_wind using "Robustness_Output.xlsx", sheet("IV_Estimates") firstrow(varlabels) cell(A1) sheetmodify

** Estimate average change in price by hour-of-day WITHOUT INSTRUMENTING but over same time period

use Setting_Sun_dataset.dta, clear
tab month, gen(m)
gen week = week(julian)
gen year_week = year*100 + week
forvalues h = 1/24 {
	quietly: reg rtm_total_lmp daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h' & solar_potential != . & wind_potential != ., vce(cluster year_week)
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..2]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..2, 1..2]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind

export excel hour-stddev_wind using "Robustness_Output.xlsx", sheet("OLS_Estimates") firstrow(varlabels) cell(A1) sheetmodify


*****************************************************************************************
** Variants of Equation (1):
** Exploiting short-run solar fluctuations (Appendix A and Figure A9)

use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab year, gen(y)
tab month, gen(m)

forvalues y = 1/4 {
	forvalues m = 1/12 {
		gen y`y'_m`m' = y`y'*m`m'
	}
}
forvalues m = 1/5 {
	gen y5_m`m' = y5*m`m'
}

forvalues h = 1/24 {
	newey rtm_total_lmp daily_solar daily_wind spot load inches_lag y1_m1-y5_m4 if hour == `h', lag(7) force
	}

	
** Day-Ahead Market Impacts (Figure A13)
	
forvalues h = 1/24 {
	quietly: newey dam_total_lmp daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..3]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..3, 1..3]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename beta_lmp3 beta_spot
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind
rename stddev_lmp3 stddev_spot

export excel hour beta_solar stddev_solar using "Robustness_Output.xlsx", sheet("DAM_Price_Change") firstrow(varlabels) cell(A1) sheetmodify

	
	
** RTM energy component impacts (Figure A4)

use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)	

forvalues h = 1/24 {
	quietly: newey rtm_energy_lmp daily_solar daily_wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..3]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..3, 1..3]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename beta_lmp3 beta_spot
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind
rename stddev_lmp3 stddev_spot

export excel hour beta_solar stddev_solar using "Robustness_Output.xlsx", sheet("RTM_Energy_Component") firstrow(varlabels) cell(A1) sheetmodify




** The impact of hourly -- instead of daily -- wind generaiton (Figure A14)

use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)

forvalues h = 1/24 {
	newey rtm_total_lmp daily_solar wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	}
	
use Setting_Sun_dataset.dta, clear
sort hour
gen id = _n
tsset id
tab month, gen(m)	

forvalues h = 1/24 {
	quietly: newey rtm_energy_lmp daily_solar wind spot load inches_lag m1-m11 if hour == `h', lag(7) force
	matrix define beta_raw = e(b)
	matrix define beta_lmp_`h' = beta_raw[1,1..3]
	matrix define cov_raw = e(V)
	matrix define cov_raw = cov_raw[1..3, 1..3]
	matrix define std_raw = cholesky(cov_raw)
	matrix define stddev_lmp_`h' = vecdiag(std_raw)
}

matrix define beta_lmp = (beta_lmp_1)
matrix define stddev_lmp = (stddev_lmp_1)

forvalues h = 2/24 {
	matrix define beta_lmp = (beta_lmp \ beta_lmp_`h')
	matrix define stddev_lmp = (stddev_lmp \ stddev_lmp_`h')
}

clear
set obs 24
gen hour = _n
svmat beta_lmp
svmat stddev_lmp
rename beta_lmp1 beta_solar
rename beta_lmp2 beta_wind
rename beta_lmp3 beta_spot
rename stddev_lmp1 stddev_solar
rename stddev_lmp2 stddev_wind
rename stddev_lmp3 stddev_spot

export excel hour beta_wind stddev_wind using "Robustness_Output.xlsx", sheet("Hourly_Wind_Impact") firstrow(varlabels) cell(A1) sheetmodify


** Nonlinear response to short-run variation in solar generation

use Setting_Sun_dataset.dta, clear
keep if year == 2016
keep if hour == 12
tab month, gen(m)

reg daily_solar daily_wind load spot inches_lag m1-m12
predict daily_solar_resid, resid

reg rtm_total_lmp daily_wind load spot inches_lag m1-m12
predict lmp_resid, resid

** Figure A10
histogram daily_solar_resid, width(2)	

sum daily_solar_resid, detail
scalar define min = r(p5)
scalar define max = r(p95)

** Figure A11
lpoly lmp_resid daily_solar_resid if daily_solar_resid>=min & daily_solar_resid<=max, noscatter bwidth(5) generate(knot fitted)


*******************************************************************************************************	
** The following section predicts counterfactual 2016 DAM prices for different levels of installed solar
** capacity. Using the predicted counterfactual prices, which are displayed in Figure A20, we then predict the operating profits for
** hypothetical conventional generators with a range of different marginal generation costs (Figure A21).

use "seasonal_dataset_temp.dta", clear

** Create daily wind and solar potential capacity factor measures
gen solar_cap_factor = solar_potential / (24 * solar_cap/1000)
gen wind_cap_factor = wind_potential / (24 * wind_cap/1000)

** Create hourly interactions for subsequent DAM price model
forvalues h = 1/24 {
	gen daily_solar_h`h' = 0
	replace daily_solar_h`h' = daily_solar if hour == `h'
}

forvalues h = 1/24 {
	gen daily_wind_h`h' = 0
	replace daily_wind_h`h' = daily_wind if hour == `h'
}

forvalues h = 1/24 {
	gen load_h`h' = 0
	replace load_h`h' = load if hour == `h'
}

forvalues h = 1/24 {
	gen spot_h`h' = 0
	replace spot_h`h' = spot if hour == `h'
}

forvalues h = 1/24 {
	gen inches_lag_h`h' = 0
	replace inches_lag_h`h' = inches_lag if hour == `h'
}

forvalues m = 1/12 {
	forvalues h = 1/24 {
		gen m`m'_h`h' = m`m' * h`h'
	}
}

** Create countefactual 2016 datasets with 2, 6, & 10 GWs of installed solar capacity

save "full_temp.dta", replace
gen solar_sample = "actual"

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 2000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_2k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 6000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_6k" if solar_sample == ""

append using "full_temp.dta"
replace solar_potential = solar_cap_factor * 10000 * (24/1000) if solar_sample == ""
replace solar_sample = "solar_10k" if solar_sample == ""

** Estimate Curtailed Solar as a function of the daily load, solar potential, wind, etc.

gen solar_load = solar_potential * daily_load
tobit solar_curtailed solar_potential solar_load daily_wind spot daily_load inches_lag m1-m11 if hour == 12 & solar_sample == "actual", ll(0)
predict solar_curtailed_fitted, ystar(0,.)

** Counterfactual daily solar is equal to the potential output minus predicted curtailment
replace daily_solar = solar_potential - solar_curtailed_fitted if solar_sample != "actual"

forvalues h = 1/24 {
	replace daily_solar_h`h' = daily_solar if hour == `h' & solar_sample != "actual"
}

** Predict counterfactual DAM prices

gen dam_fitted = .
gen resid = .

forvalues s = 1/4 {
	quietly: reg dam_total_lmp daily_solar_h1-daily_solar_h24 daily_wind_h1-daily_wind_h24 spot_h1-spot_h24 load_h1-load_h24 inches_lag_h1 inches_lag_h1-inches_lag_h24 h1-h24 if solar_sample == "actual" & season == `s'
	predict dam_fitted_`s'
	predict resid_`s', resid
	replace dam_fitted = dam_fitted_`s' if season == `s'
	replace resid = resid_`s' if solar_sample == "actual" & season == `s'
}

replace resid = resid[_n-38685] if solar_sample == "solar_2k"
replace resid = resid[_n-38685] if solar_sample == "solar_6k"
replace resid = resid[_n-38685] if solar_sample == "solar_10k"

gen dam_hat = dam_fitted + resid
drop if year != 2016
save "counterfactual_price_by_solar_temp.dta", replace

** Generate 2016 average hourly price profile by solar capacity

reg dam_hat h1-h24 if solar_sample == "actual", noc
matrix define beta_actual = e(b)'
reg dam_hat h1-h24 if solar_sample == "solar_2k", noc
matrix define beta_2k = e(b)'
reg dam_hat h1-h24 if solar_sample == "solar_6k", noc
matrix define beta_6k = e(b)'
reg dam_hat h1-h24 if solar_sample == "solar_10k", noc
matrix define beta_10k = e(b)'

matrix define avg_price = (beta_actual, beta_2k, beta_6k, beta_10k)
clear
set obs 24
gen hour = _n
svmat avg_price
rename avg_price1 actual
rename avg_price2 solar_2k
rename avg_price3 solar_6k
rename avg_price4 solar_10k

export excel hour-solar_10k using "DAM_Counterfactual.xlsx", sheet("Price_Solar") firstrow(varlabels) sheetmodify cell(A1)

** Generate kernel-smoothed distribution of predicted hourly prices

use "counterfactual_price_by_solar_temp.dta", clear
gen p = _n - 6
replace p = . if p > 110

quietly: kdensity dam_hat if solar_sample == "actual", at(p) gen(dist_actual)
quietly: kdensity dam_hat if solar_sample == "solar_2k", at(p) gen(dist_2k)
quietly: kdensity dam_hat if solar_sample == "solar_6k", at(p) gen(dist_6k)
quietly: kdensity dam_hat if solar_sample == "solar_10k", at(p) gen(dist_10k)

keep p dist_actual dist_2k dist_6k dist_10k
drop if p == .

export excel p-dist_10k using "DAM_Counterfactual.xlsx", sheet("Distribution_Solar") firstrow(varlabels) sheetmodify cell(A1)

** Predict share of hours operating and operating profits for hypothetical technologies with different marginal costs.

use "counterfactual_price_by_solar_temp.dta", clear

forvalues i = 0/100 {

	gen mc_`i' = 0
	replace mc_`i' = 1 if dam_hat >= `i'
	gen pi_`i' = 0
	replace pi_`i' = dam_hat - `i' if mc_`i' == 1
	
}

order mc_0-pi_100, sequential
collapse (sum) mc_0-pi_100, by(solar_sample)

mkmat mc_0-mc_100, matrix(operate)
mkmat pi_0-pi_100, matrix(profit)
matrix define profit = profit'
matrix define operate = operate'

clear
set obs 101
gen mc = _n - 1
svmat operate
svmat profit

foreach i in operate profit {
	rename `i'1 `i'_actual
	rename `i'2 `i'_10k
	rename `i'3 `i'_2k
	rename `i'4 `i'_6k
}

export excel mc-profit_6k using "DAM_Counterfactual.xlsx", sheet("Profit_Solar") firstrow(varlabels) sheetmodify cell(A1)
